SUBROUTINE RDSUPR(NCID)

USE PARKIND1  ,ONLY : JPIM     ,JPRB,  JPRD
USE YOMHOOK   ,ONLY : LHOOK    ,DR_HOOK, JPHOOK
USE YOMGP1S0 , ONLY : GP0      ,TSLNU0   ,TILNU0   ,QLINU0 &
                     &,FSNNU0   ,TSNNU0   ,ASNNU0   ,RSNNU0,WSNNU0 &
                     &,TRENU0   ,WRENU0 &
                     &,TLICENU0,TLMNWNU0,TLWMLNU0,TLBOTNU0,TLSFNU0 & ! FLAKE
                     &,HLICENU0,HLMLNU0 &                             ! FLAKE 
                     &,LAINU0 , BSTRNU0, BSTR2NU0,UONU0,VONU0,TONU0,SONU0 

USE YOMDPHY  , ONLY : NCSS,NLON, NLAT, NPOI, NGPP, NLALO,NVHILO,NCSNEC,NPOIP,NPOIPALL,NPOIOFF,NPOIALL

USE YOETHF   , ONLY : RHOH2O
USE YOMGC1S  , ONLY : LMASK
USE YOESOIL1S, ONLY : RDAW
USE YOMLUN1S , ONLY : NULOUT
USE YOMLOG1S , ONLY : NDIMCDF
USE YOEPHY   , ONLY : LECTESSEL, LECOLDSTART, NCWS
USE NETCDF
USE NETCDF_UTILS, ONLY: NCERROR
USE MPL_MODULE
#ifdef DOC
! (C) Copyright 2000- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**** *RDSUPR * - Reading netCDF file containing surface prognostic fields

!     Purpose.
!     --------
!            initialization surface prognostics

!**   Interface.
!     ----------
!        *CALL* *RDSUPR(NCID)

!        Explicit arguments :
!        --------------------
!     NCID      INT     NetCDF file ID


!        Implicit arguments :
!        --------------------


!     Method.
!     -------
!       Opens a file called 'soilinit' to read relevant fields
!       In the file, fields are assumed to be stored as 
!           FIELD(LAT,LON)
!       or as
!           FIELD(NLEVS,LAT,LON)


!     Externals.
!     ----------
!       NETCDF-utilities

!     Reference.
!     ----------

!     Author.
!     -------
!        Bart vd Hurk, KNMI

!     Modifications.
!     --------------
!        Original : 2000-07-13
!        E. Dutra 07/2014 : netcdf4 

!     ------------------------------------------------------------------

#endif
IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN) :: NCID



!* LOCAL VARIABLES

INTEGER ISTART2(2),ICOUNT2(2)
INTEGER ISTART3(3),ICOUNT3(3)
INTEGER NDIMS(100)
REAL(KIND=JPRB),ALLOCATABLE :: ZREALD(:)
REAL(KIND=JPRB),ALLOCATABLE :: ZBUF(:)
REAL(KIND=JPRB),ALLOCATABLE :: ZREAL3D(:,:)
REAL(KIND=JPRB),ALLOCATABLE :: ZREAL3D_DUMMY(:,:)
INTEGER NIDLON,NIDLAT,IERR,NILON,NILAT,NIDLEVS,NILEVS,NVARID,NVARTYP, &
    &    NVDIMS,NVATTS,NIDVTYPES,NIVTYPES,JK,NISNLEVS
INTEGER(KIND=JPIM) :: NMX,NMY,ILEVS,JL,J,IVTYPES,JD,DIMLEN,NIMON
CHARACTER*100 CNAME
CHARACTER*20 CDUM

LOGICAL LOPRES,LRSNML,LEXTENDSOIL
INTEGER IVAR,NDIM,NGATTS,IRECDIM,NVARS,STATUS
INTEGER(KIND=JPIM)           :: NVARS2D,NVARS3D
CHARACTER(LEN=20)            :: CVARS2D(40),CVARS3D(40),CVAR
INTEGER(KIND=JPIM)           :: MYPROC, NPROC
INTEGER(KIND=JPIM)           :: ISP, IENP
REAL(KIND=JPHOOK)              :: ZHOOK_HANDLE

DATA ISTART2/1,1/
DATA ISTART3/1,1,1/

#include "minmax.intfb.h"

IF (LHOOK) CALL DR_HOOK('RDSUPR',0,ZHOOK_HANDLE)

MYPROC = MPL_MYRANK()
NPROC  = MPL_NPROC()

ISP   = NPOIOFF(MYPROC)+1
IENP = NPOIOFF(MYPROC+1)

NILON=-999
NILAT=-999
NILEVS=-999
NIVTYPES=-999
NISNLEVS=-999


IF( MYPROC == 1 ) THEN
!* check the dimensions
  CALL NCERROR (NF90_INQUIRE(NCID, nDimensions=NDIM,nVariables=NVARS),'FINDING N DIMS')

  DO JD=1,NDIM
    CALL NCERROR(NF90_INQUIRE_DIMENSION(NCID, JD, name=CNAME, len=DIMLEN),'INQ DIMENSION')
    SELECT CASE (TRIM(CNAME))
      CASE('x','lon')
        NILON=DIMLEN
        WRITE(NULOUT,*)'DIMENSION x ',TRIM(CNAME),' FOUND WITH LEN=',NILON
     CASE('y','lat')
        NILAT=DIMLEN
        WRITE(NULOUT,*)'DIMENSION y ',TRIM(CNAME),' FOUND WITH LEN=',NILAT
     CASE('nlevs')
        NILEVS=DIMLEN
        WRITE(NULOUT,*)'DIMENSION nlevs ',TRIM(CNAME),' FOUND WITH LEN=',NILEVS
     CASE('vtype')
        NIVTYPES=DIMLEN
        WRITE(NULOUT,*)'DIMENSION vtype ',TRIM(CNAME),' FOUND WITH LEN=',NIVTYPES
     CASE('nlevsn')
        NISNLEVS=DIMLEN
        WRITE(NULOUT,*)'DIMENSION vtype ',TRIM(CNAME),' FOUND WITH LEN=',NISNLEVS
    END SELECT
  END DO
ENDIF

CALL MPL_BROADCAST(NDIM,    KROOT=1,KTAG=100,CDSTRING='NDIM')
CALL MPL_BROADCAST(NILON,   KROOT=1,KTAG=200,CDSTRING='NILON')
CALL MPL_BROADCAST(NILAT,   KROOT=1,KTAG=300,CDSTRING='NILAT')
CALL MPL_BROADCAST(NILEVS,  KROOT=1,KTAG=400,CDSTRING='NILEVS')
CALL MPL_BROADCAST(NIVTYPES,KROOT=1,KTAG=500,CDSTRING='NIVTYPES')
CALL MPL_BROADCAST(NISNLEVS,KROOT=1,KTAG=600,CDSTRING='NISNLEVS')

IF (NDIMCDF == 1 ) THEN
  NILAT = 0
ENDIF
LRSNML=.TRUE.
LEXTENDSOIL=.FALSE.
IF(NILON /= NLON .OR. NILAT /= NLAT) THEN
  WRITE(NULOUT,*)'NLON OR NLAT NOT SPECIFIED CORRECTLY:'
  WRITE(NULOUT,*)'NLON IN NAMELIST AND SOILINIT: ',NLON,NILON
  WRITE(NULOUT,*)'NLAT IN NAMELIST AND SOILINIT: ',NLAT,NILAT
  CALL ABOR1('RDSUPR:')
ENDIF
IF(NILEVS /= NCSS) THEN
  WRITE(NULOUT,*)'NLEVS NOT SPECIFIED CORRECTLY:'
  WRITE(NULOUT,*)'NLEVS IN CODE AND SOILINIT: ',NCSS,NILEVS
  IF(NILEVS == 4 .AND. NCSS > 4) THEN
     WRITE(NULOUT,*)'Starting from the 4 layer scheme'
     WRITE(NULOUT,*)'All levels below layer 4 set equal to layer 4, this should be just a temporary fix'
     LEXTENDSOIL=.TRUE.
  ELSE
    CALL ABOR1('RDSUPR:')
  ENDIF
ENDIF
IF (NIVTYPES == -999 ) THEN
  NIVTYPES=NVHILO
  WRITE(NULOUT,*)'vtype DIMENSION NOT FOUND IN NETCDF FILE!!! SET==NVHILO'
ENDIF
IF(NIVTYPES /= NVHILO) THEN
  WRITE(NULOUT,*)'NVHILO NOT SPECIFIED CORRECTLY:'
  WRITE(NULOUT,*)'NVHILO IN CODE AND PROGINIT: ',NVHILO,NIVTYPES
  CALL ABOR1('RDSUPR:')
ENDIF
IF (NISNLEVS == -999  ) THEN
  IF  (NCSNEC>1) THEN
    WRITE(NULOUT,*)'Input does not contain snow multi-layer'
    WRITE(NULOUT,*) 'Input initialized from total fields'
    LRSNML=.FALSE.
  ELSE
    LRSNML=.FALSE.
  ENDIF
ENDIF
IF(NISNLEVS /= -999 .AND. (NCSNEC /= NISNLEVS) ) THEN
  WRITE(NULOUT,*)'Input contains:',NISNLEVS,' snow layer, while running with:',NCSNEC,' layers'
  WRITE(NULOUT,*) 'Input initialized from total fields'
  LRSNML=.FALSE.
ENDIF




!* Note the reversal of the order of array storage

IF(NDIMCDF == 2)THEN
  ICOUNT2(1)=NILON
  ICOUNT2(2)=NILAT

  ICOUNT3(1)=NILON
  ICOUNT3(2)=NILAT
  ICOUNT3(3)=NILEVS
  NMX=NLON
  NMY=NLAT
ELSE
  ICOUNT2(1)=NILON

  ICOUNT3(1)=NILON
  ICOUNT3(2)=NILEVS
  NMX=NLALO
  NMY=1
ENDIF


!* Read surface prognostic values

! 3D FIELDS
ALLOCATE (ZREAL3D(NLALO,NCSS))
ALLOCATE (ZREAL3D_DUMMY(NLALO,NCSS))
ALLOCATE (ZBUF(NPOIALL))
NVARS3D=3
 CVARS3D(1:NVARS3D)=(/'SoilTemp ','icetemp  ','SoilMoist'   /)
DO IVAR=1,NVARS3D
  CVAR=TRIM(CVARS3D(IVAR))
  IF( MYPROC == 1 ) THEN
    STATUS = NF90_INQ_VARID(NCID, CVAR, NVARID)
    IF ( STATUS /= 0 ) THEN
      WRITE(NULOUT,'(A)') CVAR//' NOT PRESENT IN INPUT FILE'
    ELSE
      CALL NCERROR( NF90_GET_VAR(NCID,NVARID,ZREAL3D,ISTART3,ICOUNT3),'READING '//CVAR)
    ENDIF
    IF (LEXTENDSOIL) THEN
    ZREAL3D_DUMMY=ZREAL3D
     DO ILEVS=1,3
       ZREAL3D(:,ILEVS) = ZREAL3D_DUMMY(:,1)
     ENDDO
     DO ILEVS=4,5
       ZREAL3D(:,ILEVS) = ZREAL3D_DUMMY(:,2)
     ENDDO
     DO ILEVS=6,7
       ZREAL3D(:,ILEVS) = ZREAL3D_DUMMY(:,3)
     ENDDO
     DO ILEVS=8,NCSS
       ZREAL3D(:,ILEVS) = ZREAL3D_DUMMY(:,4)
     ENDDO
     DEALLOCATE(ZREAL3D_DUMMY)
    ENDIF
  ENDIF
  DO ILEVS=1,NCSS
    SELECT CASE(CVAR)
      CASE('SoilTemp')
        IF( MYPROC == 1 ) THEN
          IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
        ENDIF
        CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREAL3D(:,ILEVS),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:TSLNU0")
        TSLNU0(:,ILEVS)=PACK(ZBUF(:),LMASK(ISP:IENP))
      CASE('SoilMoist')
        IF( MYPROC == 1 ) THEN
          IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
          CALL NCERROR( NF90_GET_ATT(NCID, NVARID, 'units', CNAME),'Getting SM units')
          IF(CNAME(1:6) /= 'm3 m-3'.AND. CNAME(1:5) /= 'm3/m3' .AND. CNAME(1:10) /= 'm**3 m**-3') THEN
            ZREAL3D(:,ILEVS)=ZREAL3D(:,ILEVS)/(RHOH2O*RDAW(ILEVS))
            WRITE(NULOUT,*)'RDSUPR: INPUT SOIL MOISTURE UNITS:',CNAME(1:10),' CONVERTED TO M3 M-3'
          ENDIF
        ENDIF
        CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREAL3D(:,ILEVS),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:QLINU0")
        QLINU0(:,ILEVS)=PACK(ZBUF(:),LMASK(ISP:IENP))
      CASE('icetemp')
        IF( MYPROC == 1 ) THEN
          IF ( STATUS /= 0 ) ZREAL3D(:,ILEVS) = 273._JPRB
        ENDIF
        CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREAL3D(:,ILEVS),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:TILNU0")
        TILNU0(:,ILEVS)=PACK(ZBUF(:),LMASK(ISP:IENP))
      CASE DEFAULT
        WRITE(NULOUT,*) CVAR, ' Not defined in RDCLIM'
        CALL ABOR1('RDSUPR:')
    END SELECT
    IF( MYPROC == 1 ) THEN
      WRITE(CDUM,'(A8,I2.2)')TRIM(CVAR),ILEVS
      CALL MINMAX(CDUM,ZREAL3D(:,ILEVS),NMX,NMY,LMASK,NULOUT)
    ENDIF
    CALL MPL_BARRIER()
  ENDDO
ENDDO


ALLOCATE (ZREALD(NLALO))
NVARS2D=14
 CVARS2D(1:NVARS2D)=(/'SWE      ','snowdens ','SAlbedo  ','SnowT    ',&
                      'CanopInt ','AvgSurfT ','TLICE    ','TLMNW    ',&
                      'TLWML    ','TLBOT    ','TLSF     ','HLICE    ',&
                      'HLML     ','slw      '/)

DO IVAR=1,NVARS2D
  CVAR=TRIM(CVARS2D(IVAR))
  IF( MYPROC == 1 ) THEN
    STATUS = NF90_INQ_VARID(NCID, CVAR, NVARID)
    IF ( STATUS /= 0 ) THEN
      WRITE(NULOUT,'(A)') CVAR//' NOT PRESENT IN INPUT FILE'
    ELSE
      CALL NCERROR( NF90_GET_VAR(NCID,NVARID,ZREALD,ISTART2,ICOUNT2),'READING '//CVAR)
    ENDIF
  ENDIF
  SELECT CASE(CVAR)
    CASE('SWE')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:FSNNU0")
      FSNNU0(:,1)=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('snowdens')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:RSNNU0")
      RSNNU0(:,1)=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('SAlbedo')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:ASNNU0")
      ASNNU0=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('SnowT')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:TSNNU0")
      TSNNU0(:,1)=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('CanopInt')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:WRENU0")
      WRENU0=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('AvgSurfT')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:TRENU0")
      TRENU0=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('TLICE')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) ZREALD(:) = 273._JPRB
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:TLICENU0")
      TLICENU0=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('TLMNW')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) ZREALD(:) = 288._JPRB 
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:TLMNWNU0")
      TLMNWNU0=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('TLWML')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) ZREALD(:) = 288._JPRB
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:TLWMLNU0")
      TLWMLNU0=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('TLBOT')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) ZREALD(:) = 288._JPRB
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:TLBOTNU0")
      TLBOTNU0=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('TLSF')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) ZREALD(:) = 0.65_JPRB
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:TLSFNU0")
      TLSFNU0=PACK(ZBUF(:),LMASK(ISP:IENP)) 
    CASE('HLICE')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) ZREALD(:) = 0._JPRB
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:HLICENU0")
      HLICENU0=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('HLML')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) ZREALD(:) = 10._JPRB
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:HLMLNU0")
      HLMLNU0=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE('slw')
      IF( MYPROC == 1 ) THEN
        IF ( STATUS /= 0 ) ZREALD(:) = 0._JPRB
        WRITE(NULOUT,'(A)') CVAR//' Set to Zero'
      ENDIF
      CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREALD(:),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:WSNNU0")
      WSNNU0(:,1)=PACK(ZBUF(:),LMASK(ISP:IENP))
    CASE DEFAULT
      WRITE(NULOUT,*) CVAR, ' Not defined in RDCLIM'
      CALL ABOR1('RDSUPR:')
  END SELECT
  IF( MYPROC == 1 ) THEN
    WRITE(CDUM,*)TRIM(CVAR)
    CALL MINMAX(CDUM,ZREALD,NMX,NMY,LMASK,NULOUT)
  ENDIF
  CALL MPL_BARRIER()
ENDDO

!! Read snow ML fields
IF (( LRSNML ) .AND. (.NOT. LECOLDSTART) ) THEN
  IF(NDIMCDF == 2)THEN
    ICOUNT2(1)=NILON
    ICOUNT2(2)=NILAT

    ICOUNT3(1)=NILON
    ICOUNT3(2)=NILAT
    ICOUNT3(3)=NCSNEC
    NMX=NLON
    NMY=NLAT
  ELSE
    ICOUNT2(1)=NILON

    ICOUNT3(1)=NILON
    ICOUNT3(2)=NCSNEC
    NMX=NLALO
    NMY=1
  ENDIF
  DEALLOCATE (ZREAL3D)
  ALLOCATE (ZREAL3D(NLALO,NCSNEC))
  
  NVARS3D=4
  CVARS3D(1:NVARS3D)=(/'SWEML     ','snowdensML','SnowTML   ','slwML     '/)
  DO IVAR=1,NVARS3D
    CVAR=TRIM(CVARS3D(IVAR))
    IF( MYPROC == 1 ) THEN
      STATUS = NF90_INQ_VARID(NCID, CVAR, NVARID)
      IF ( STATUS /= 0 ) THEN
        WRITE(NULOUT,'(A)') CVAR//' NOT PRESENT IN INPUT FILE'
      ELSE
        CALL NCERROR( NF90_GET_VAR(NCID,NVARID,ZREAL3D,ISTART3,ICOUNT3),'READING '//CVAR)
      ENDIF
    ENDIF
    DO ILEVS=1,NCSNEC
      SELECT CASE(CVAR)
        CASE('SWEML')
          IF( MYPROC == 1 ) THEN
            IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
          ENDIF
          CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREAL3D(:,ILEVS),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:FSNNU0")
          FSNNU0(:,ILEVS)=PACK(ZBUF(:),LMASK(ISP:IENP))
        CASE('snowdensML')
          IF( MYPROC == 1 ) THEN
            IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
          ENDIF
          CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREAL3D(:,ILEVS),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:RSNNU0")
          RSNNU0(:,ILEVS)=PACK(ZBUF(:),LMASK(ISP:IENP))
        CASE('SnowTML')
          IF( MYPROC == 1 ) THEN
            IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
          ENDIF
          CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREAL3D(:,ILEVS),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:TSNNU0")
          TSNNU0(:,ILEVS)=PACK(ZBUF(:),LMASK(ISP:IENP))
        CASE('slwML')
          IF( MYPROC == 1 ) THEN
            IF ( STATUS /= 0 ) CALL ABOR1('RDSUPR:')
          ENDIF
          CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREAL3D(:,ILEVS),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:WSNNU0")
          WSNNU0(:,ILEVS)=PACK(ZBUF(:),LMASK(ISP:IENP))
        CASE DEFAULT
          WRITE(NULOUT,*) CVAR, ' Not defined in RDCLIM'
          CALL ABOR1('RDSUPR:')
      END SELECT
      IF( MYPROC == 1 ) THEN
        WRITE(CDUM,'(A8,I2.2)')TRIM(CVAR),ILEVS
        CALL MINMAX(CDUM,ZREAL3D(:,ILEVS),NMX,NMY,LMASK,NULOUT)
      ENDIF
      CALL MPL_BARRIER()
    ENDDO
  ENDDO
  
ELSE
  ! Fill up snow 3D fields 
  WRITE(NULOUT,'(A)') 'Filling Snow 3D fields with 2D information'
  DO JK=2,NCSNEC
    FSNNU0(:,JK) = 0.0_JPRB
    WSNNU0(:,JK) = 0.0_JPRB
    RSNNU0(:,JK) = RSNNU0(:,1)
    TSNNU0(:,JK) = TSNNU0(:,1)
  ENDDO

ENDIF 
    



!* < CTESSEL
!    IF (LECTESSEL) THEN 
DEALLOCATE (ZREAL3D)
ALLOCATE (ZREAL3D(NLALO,NVHILO))

IF(NDIMCDF == 2)THEN
  ICOUNT3(1)=NILON
  ICOUNT3(2)=NILAT
  ICOUNT3(3)=NIVTYPES
ELSE
  ICOUNT3(1)=NILON
  ICOUNT3(2)=NIVTYPES
ENDIF

NVARS3D=3
 CVARS3D(1:NVARS3D)=(/'laivt    ','Biomstr  ','Biomstr2 '   /)
DO IVAR=1,NVARS3D
  CVAR=TRIM(CVARS3D(IVAR))
  IF( MYPROC == 1 ) THEN
    STATUS = NF90_INQ_VARID(NCID, CVAR, NVARID)
    IF ( STATUS /= 0 ) THEN
      WRITE(NULOUT,'(A)') CVAR//' NOT PRESENT IN INPUT FILE'
    ELSE
      CALL NCERROR( NF90_GET_VAR(NCID,NVARID,ZREAL3D,ISTART3,ICOUNT3),'READING '//CVAR)
    ENDIF
  ENDIF
  DO IVTYPES=1,NVHILO
    SELECT CASE(CVAR)
      CASE('laivt')
        IF( MYPROC == 1 ) THEN
          IF ( STATUS /= 0 ) ZREAL3D(:,IVTYPES) = 0.3_JPRB
        ENDIF
        CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREAL3D(:,IVTYPES),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:LAINU0")
        LAINU0(:,IVTYPES)=PACK(ZBUF(:),LMASK(ISP:IENP))
      CASE('Biomstr')
        IF( MYPROC == 1 ) THEN
          IF ( STATUS /= 0 ) ZREAL3D(:,IVTYPES) = 0._JPRB
        ENDIF
        CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREAL3D(:,IVTYPES),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:BSTRNU0")
        BSTRNU0(:,IVTYPES)=PACK(ZBUF(:),LMASK(ISP:IENP))
      CASE('Biomstr2')
        IF( MYPROC == 1 ) THEN
          IF ( STATUS /= 0 ) ZREAL3D(:,IVTYPES) = 0._JPRB
        ENDIF
        CALL MPL_SCATTERV(PRECVBUF=ZBUF(:),KROOT=1,PSENDBUF=ZREAL3D(:,IVTYPES),KSENDCOUNTS=NPOIPALL(:),CDSTRING="RDSUPR:BSTR2NU0")
        BSTR2NU0(:,IVTYPES)=PACK(ZBUF(:),LMASK(ISP:IENP))
      CASE DEFAULT
        WRITE(NULOUT,*) CVAR, ' Not defined in RDCLIM'
        CALL ABOR1('RDSUPR:')
    END SELECT
    IF( MYPROC == 1 ) THEN
      WRITE(CDUM,'(A8,I2.2)')TRIM(CVAR),ILEVS
      CALL MINMAX(CDUM,ZREAL3D(:,IVTYPES),NMX,NMY,LMASK,NULOUT)
    ENDIF
    CALL MPL_BARRIER()
  ENDDO
ENDDO

!* Set ocean data to zero - this should be removed at some point ! 
UONU0(:,:) = 0._JPRB
VONU0(:,:) = 0._JPRB
TONU0(:,:) = 0._JPRB
SONU0(:,:) = 0._JPRB

!* Check final values

DO JL=1,NPOI
  DO J=1,NGPP
    IF (GP0(JL,J) == -999._JPRB) THEN
      WRITE(NULOUT,*) 'ERROR SOIL INITIALISATION', JL,J
      CALL ABOR1('RDSUPR:')
    ENDIF
  ENDDO
ENDDO
DEALLOCATE (ZREAL3D)
DEALLOCATE (ZREALD)

CALL MPL_BARRIER()

IF (LHOOK) CALL DR_HOOK('RDSUPR',1,ZHOOK_HANDLE)

END SUBROUTINE RDSUPR
